A previous presentation considered an analytic approximation to the planar orbits generated by a Hamiltonian of the form

$H=\frac{p^2}{2m}+\alpha r^k=E$

That approximation made heavy use of Jacobi elliptic functions and elliptic integrals. While these are standard special functions appearing regularly in mathematical physics, the overall approximation is perhaps not as clear as one would like. There is indeed a much simpler approximation whose parameters can be related simply to the behavior of exact orbits.

The equations of motion in Cartesian coordinates are

$\begin{array}{l}\stackrel{··}{x}=\frac{{\stackrel{·}{p}}_x}{m}=-\frac{1}{m}\frac{\partial H}{\partial x}=-\frac{k\alpha }{m}{r}^{k-2}x\\ \stackrel{··}{y}=\frac{{\stackrel{·}{p}}_y}{m}=-\frac{1}{m}\frac{\partial H}{\partial y}=-\frac{k\alpha }{m}{r}^{k-2}y\end{array}$

Taking $m=\alpha =1$ for convenience and starting the orbit from a turning point, the parametric behavior for a variety of inputs looks like this:

For $n=2$ the orbit is a closed ellipse for all other inputs. That immediately raises the question as to whether the orbits for other powers are approximately elliptic as well, but rotating. That is, start with an ellipse of the form

$x\left(t\right)=acos\omega ty\left(t\right)=bsin\omega t$

and give it a temporally dependent rotation:

$\left[\begin{array}{c}X\left(t\right)\\ Y\left(t\right)\end{array}\right]=\left[\begin{array}{cc}cos\Delta t& -sin\Delta t\\ sin\Delta t& cos\Delta t\end{array}\right]\left[\begin{array}{c}x\left(t\right)\\ y\left(t\right)\end{array}\right]=\left[\begin{array}{c}acos\omega tcos\Delta t-bsin\omega tsin\Delta t\\ acos\omega tsin\Delta t+bsin\omega tcos\Delta t\end{array}\right]$

Under this rotation, the radial combination of the two Cartesian coordinates is naturally unchanged:

$X^2\left(t\right)+Y^2\left(t\right)=a^2{cos}^2\omega t+b^2{sin}^2\omega t=x^2\left(t\right)+y^2\left(t\right)$

Using trigonometric product rules, the individual Cartesian coordinates can be written

$\begin{array}{l}X\left(t\right)=\frac{a+b}{2}cos(\omega +\Delta )t+\frac{a-b}{2}cos(\omega -\Delta )t\\ Y\left(t\right)=\frac{a+b}{2}sin(\omega +\Delta )t-\frac{a-b}{2}sin(\omega -\Delta )t\end{array}$

These rotated ellipses are no longer exactly centered on the origin, which can be seen by setting $\Delta =\omega$ in the first coordinate and noting the constant offset along the x-axis.

The choice of one parameter is immediately clear: the orbits must start at the same point. That means $a=r_0$ identically.

This parametrization of the orbit can be explored with this interactive graphic:

Manipulating the last two inputs rotates the parametrized curve in either direction. This makes sense because the parametrization depends upon linear combinations of these two parameters, so they should work somewhat similarly. This is also somewhat misleading because there appears to be multiple ways to make the parametrization fit the actual orbits. This is a figment of this particular visualization, since there are two different times involved in the orbits, and will be remedied by visualizing individual coordinates.

Manipulating the fourth input, for the parameter b, expands or flattens the rotated ellipse. This indicates the way to determine this parameter: it should be set equal to the second turning point of the orbit. For $k=2$ this second turning point is exactly the semiminor axis of the nonrotated ellipse, and for the general case it can be determined numerically. Writing the Hamiltonian in polar coordinates,

$H=\frac{p_r^2}{2m}+\frac{L^2}{2m{r}^2}+\alpha r^k=E$

one can determine the energy from the first turning point, where radial momentum goes to zero:

$E=\frac{L^2}{2m{r}_0^2}+\alpha r_0^k$

The second turning point is then the solution to the nonlinear equation

$\begin{array}{l}2mE{r}_1^2-2ma{r}_1^{k+2}-L^2=0\end{array}$

where multiplying through by the angular momentum denominator improves numerical stability of root finding. With the choice $b=r_1$ , the visualization mentioned above for the X-coordinate is

while that for the Y-coordinate is

Now manipulating the last two inputs indicates there is only one way for the parametrization to match the exact orbits. The difference in their absolute values indicates how to distinguish between them.

There are two important integrals in the context of power potential orbits, what may be called the energy integral and the angular momentum integral. One nice way to derive these integrals is through Hamilton-Jacobi theory. The time taken in moving between the two turning points is given by the energy integral:

$\Delta t=\underset{r_1}{\overset{r_0}{\int }}\frac{mrdr}{\sqrt{2mE{r}^2-2m\alpha r^{k+2}-L^2}}$

The direction of integration produces a positive result. The period of the orbit can be taken to be four times this value. It is interesting to compare this exact period (in blue) to that of limiting circular orbits for power potentials:

Since the two are remarkably similar, the value for circular orbits has been used throughout this presentation to set the domains for integrating differential equations, as well as the physical limits on angular momentum for bound orbits.

The angular change in moving between turning points is given by the angular momentum integral:

$\Delta \phi =\underset{r_1}{\overset{r_0}{\int }}\frac{Ldr}{r\sqrt{2mE{r}^2-2m\alpha r^{k+2}-L^2}}$

For $n\ne 2$ this value differs from $\frac{\pi }{2}$ , which is why most power potential orbits are not closed. Here is the comparison:

Since the period of the orbit is a much larger number compared to the angular discrepancy, one would expect the larger of the last two parameters to correspond to exact period. From manipulating graphics above this is clearly ω, with Δ providing a smaller adjustment to the parametrization. This also makes intuitive sense: the period of the nonrotated ellipse should be the baseline time frame matching the exact orbit.

For numerical correspondence, the period of the underlying circular functions needs to be adjusted to that of the exact period, meaning

$\omega =\frac{2\pi }{4\Delta t}=\frac{\pi }{2\Delta t}$

This choice correlates with the choice of b, as can be seen from a simple evaluation:

$r_1^2=r^2\left(\Delta t\right)=a^2{cos}^2\left(\frac{\pi }{2\Delta t}\Delta t\right)+b^2{sin}^2\left(\frac{\pi }{2\Delta t}\Delta t\right)=b^2$

The choices for a and b guarantee that the approximation matches the exact result at radial extrema, and the choice for ω sets the basic period of the approximation. The final parameter must then be used to guarantee the approximation matches the exact result for angular changes, and as such must be proportional to the angular discrepancy. It can be evaluated as

$\Delta =\omega (\frac{\Delta \phi }{\pi /2}-1)=\frac{1}{\Delta t}(\Delta \phi -\frac{\pi }{2})$

Including these the choices for the final two parameters, the visualization for the X-coordinate is

and that for the Y-coordinate is

In both cases the approximation improves for larger angular momentum, which is not surprising as these will be more circular orbits. The approximation is less accurate for low angular momentum, but given its simplicity that is also not surprising. In this same regime the approximation from the previous presentation is more accurate, but again more complicated mathematically.

In comparing the two approximations, it is at first curious to note that

$\omega =w_{\mathrm{radial}}\Delta =-w_{\mathrm{beats}}$

until one realizes that the quantities arise from the same two fundamental integrals and thus must necessarily be related.

For negative exponents in the potential, the coupling constant must change sign for bound orbits, the remaining graphics will take $a=-1$ . Since it was noted above that the rotated ellipses do not remain centered exactly on the origin, one can attempt to blithely apply the solution above to these new cases. Proceeding in this way, the visualization for the X-coordinate is

while the corresponding one for the Y-coordinate is

Surprisingly enough these approximations are not too bad. This means that the rotation of the ellipse originally centered on the origin moves fairly close to one centered on a focus of the new ellipse. This can be taken as an approximate analytic transformation between the two types of ellipses, those for positive and negative exponents in the potential. Unfortunately a visualization of the radial variable shows the overall inaccuracy of this approximation more clearly:

The approximation for positive exponents began with an ellipse centered on the origin, since that is the exact solution for a simple harmonic oscillator. For negative exponents one should really start with an ellipse displaced to a focus, since that is the exact solution for the Kepler potential. Such an ellipse has a parametrization

$x\left(t\right)=acos\left(\omega t\right)+cy\left(t\right)=bsin\left(\omega t\right)$

where parameters are connected by $c^2=a^2-b^2$ . Values here are related to turning points by

$\begin{array}{c}a+c=r_0\\ a-c=r_1\end{array}\to \begin{array}{c}a=\frac{r_0+r_1}{2}\\ c=\frac{r_0-r_1}{2}\end{array}$

Two other alterations are needed: a typical orbit for positive exponents requires four trips between turning points for a complete period, while with negative exponents only two are required, which means the frequency is doubled. Additionally, the reference point for angular change with negative exponents is the Kepler potential, for which the appropriate angular value is also twice as large. These two alterations mean

$\omega =\frac{\pi }{\Delta t}\Delta =\frac{1}{\Delta t}(\Delta \phi -\pi )$

With these modifications to the approximation, the visualization of the radial variable is much improved:

Unfortunately, the individual Cartesian coordinates are not as much improved. First the X-coordinate,

followed by the Y-coordinate:

Clearly orbits for negative exponents are not as close to rotating ellipses as those for positive exponents. Curious...

In summary, orbits starting from a turning point for the power-potential Hamiltonian

$H=\frac{p_r^2}{2m}+\frac{L^2}{2m{r}^2}+\alpha r^k=E=\frac{L^2}{2m{r}_0^2}+\alpha r_0^k$

can be relatively accurately parametrized by the simple analytic approximation

$\left[\begin{array}{c}X\left(t\right)\\ Y\left(t\right)\end{array}\right]=\left[\begin{array}{c}{r}_{0}cos\omega tcos\Delta t-r_{1}sin\omega tsin\Delta t\\ r_{0}cos\omega tsin\Delta t+r_{1}sin\omega tcos\Delta t\end{array}\right]$

where the second turning point is the numerical solution to the nonlinear equation

$\begin{array}{l}2mE{r}_1^2-2ma{r}_1^{k+2}-L^2=0\end{array}$

With two auxiliary quadratures

$\Delta t=\underset{r_1}{\overset{r_0}{\int }}\frac{mrdr}{\sqrt{2mE{r}^2-2m\alpha r^{k+2}-L^2}}$

$\Delta \phi =\underset{r_1}{\overset{r_0}{\int }}\frac{Ldr}{r\sqrt{2mE{r}^2-2m\alpha r^{k+2}-L^2}}$

the two final parameters are

$\omega =\frac{\pi }{2\Delta t}\Delta =\frac{1}{\Delta t}(\Delta \phi -\frac{\pi }{2})$

The approximation is most accurate with positive exponents in the potential. For negative exponents one can switch to initially displaced ellipses as described above.

Uploaded 2022.10.31 — Updated 2022.11.10 analyticphysics.com